#include <gmp.h>
#include <pthread.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>

#define NUMITT 23

//memória compartilhada pelas threads
typedef struct {

	mpf_t a1, a2, a3, b1, b2, b3, t, p;
	//a<n-1>, a<n>, a<n+1>
	//b<n-1>, b<n>, b<n+1>
	//t<n>, p<n>

} shared_mem;

//cria uma memória compartilhada já inicializada
shared_mem * create_shared_mem() {

	shared_mem *sh_mem = (shared_mem*)malloc(sizeof(shared_mem));

	mpf_t aux;
	mpf_init(aux);
	mpf_set_prec(aux,33300000);
	mpf_set_d(aux,2.0); //aux = 2
	mpf_sqrt(aux,aux); //aux = sqrt(2)
	
	mpf_init(sh_mem->a3);
	mpf_set_prec(sh_mem->a3, 33300000);
	mpf_set_d(sh_mem->a3,1.0); //a3 = 1

	mpf_init(sh_mem->a2);
	mpf_set_prec(sh_mem->a2, 33300000);
	mpf_set_d(sh_mem->a2,0.0);

	mpf_init(sh_mem->a1);
	mpf_set_prec(sh_mem->a1, 33300000);
	mpf_set_d(sh_mem->a1,0.0);

	mpf_init(sh_mem->b1);
	mpf_set_prec(sh_mem->b1, 33300000);
	mpf_set_d(sh_mem->b1,1.0);

	mpf_init(sh_mem->b3);
	mpf_set_prec(sh_mem->b3, 33300000);
	mpf_set_d(sh_mem->b3,1.0);
	mpf_div(sh_mem->b3, sh_mem->b3, aux); //b3 = 1 / sqrt(2)

	mpf_clear(aux);

	mpf_init(sh_mem->b2);
	mpf_set_prec(sh_mem->b2, 33300000);
	mpf_set_d(sh_mem->b2,1.0);
	
	mpf_init(sh_mem->t);
	mpf_set_prec(sh_mem->t, 33300000);
	mpf_set_d(sh_mem->t,1.0); //t = 1
	mpf_div_ui(sh_mem->t, sh_mem->t, 4); //t = 1/4	

	mpf_init(sh_mem->p);
	mpf_set_prec(sh_mem->p, 33300000);
	mpf_set_d(sh_mem->p,1.0);	


	return sh_mem;

}

//destrói memória compartilhada liberando seu espaço de memória
void destroy_shared_mem(shared_mem * sh_mem) {

	mpf_clear(sh_mem->a1);
	mpf_clear(sh_mem->a2);
	mpf_clear(sh_mem->a3);
	mpf_clear(sh_mem->b1);
	mpf_clear(sh_mem->b2);
	mpf_clear(sh_mem->b3);
	mpf_clear(sh_mem->t);
	free(sh_mem);

}

//configurações iniciais para cada iteração
void config_iteration(shared_mem *sh_mem) {

	mpf_set(sh_mem->a1, sh_mem->a2); //a<k-1> = a<k>
	mpf_set(sh_mem->b1, sh_mem->b2); //b<k-1> = b<k>
	mpf_set(sh_mem->a2, sh_mem->a3); //a<k> = a<k+1>	
	mpf_set(sh_mem->b2, sh_mem->b3); //b<k> = b<k+1>	

}

//calcula a<n>
void *calc_a(void *args) {

	shared_mem *sh_mem = (shared_mem*) args;

	mpf_add(sh_mem->a3, sh_mem->b2, sh_mem->a2); //a<k+1> = b<k> + a<k>
	mpf_div_ui(sh_mem->a3, sh_mem->a3, 2);  //a<k+1> = (b<k> + a<k>)/2

	pthread_exit(NULL);

}

//calcula b<n>
void *calc_b(void *args) {

	shared_mem *sh_mem = (shared_mem*) args;

	mpf_mul(sh_mem->b3, sh_mem->b2, sh_mem->a2); //b<k+1> = b<k> * a<k>
	mpf_sqrt(sh_mem->b3, sh_mem->b3); //b<k+1> = sqrt(b<k> * a<k>)

	

	pthread_exit(NULL);

}

//calcula t<n>
void *calc_t(void *args) {
	
	shared_mem *sh_mem = (shared_mem*) args;

	mpf_sub(sh_mem->a1,sh_mem->a1,sh_mem->a2); //aux = a<k> - a<k+1>
	mpf_pow_ui(sh_mem->a1, sh_mem->a1, 2); //aux = (a<k> - a<k+1>)²
	mpf_mul(sh_mem->a1, sh_mem->a1, sh_mem->p); //t<k> = p(a<k> - a<k+1>)²
	mpf_sub(sh_mem->t, sh_mem->t, sh_mem->a1);
	mpf_mul_ui(sh_mem->p,sh_mem->p,2);

	pthread_exit(NULL);

}

int main() {

	FILE *f;
	mpf_t pi;
	mpf_init(pi);
	mpf_set_prec(pi,33300000);

	pthread_t threads[2];
	shared_mem *sh_mem = create_shared_mem();

	int i = 0;

	config_iteration(sh_mem);
	//calcula a<n+1> e b<n+1>
	pthread_create(&threads[0], NULL, calc_a, sh_mem);
	pthread_create(&threads[1], NULL, calc_b, sh_mem);
	pthread_join(threads[0],NULL);
	pthread_join(threads[1],NULL);

	for ( i = 0; i < NUMITT-1; i++ ) {

		config_iteration(sh_mem);
		//calcula a<n+1> e b<n+1>
		pthread_create(&threads[0], NULL, calc_a, sh_mem);
		pthread_create(&threads[1], NULL, calc_b, sh_mem);
		//calcula t<n>
		pthread_create(&threads[2], NULL, calc_t, sh_mem);
		
		pthread_join(threads[0],NULL);
		pthread_join(threads[1],NULL);
		pthread_join(threads[2],NULL);		

	}

	config_iteration(sh_mem);
	//calcula t<n+1>
	mpf_sub(sh_mem->a1,sh_mem->a1,sh_mem->a2); //a<k> = a<k> - a<k+1>
	mpf_pow_ui(sh_mem->a1, sh_mem->a1, 2); //a<k> = (a<k> - a<k+1>)²
	mpf_mul(sh_mem->a1, sh_mem->a1, sh_mem->p); //a<k> = p(a<k> - a<k+1>)²
	mpf_sub(sh_mem->t, sh_mem->t, sh_mem->a1); //t<k> = t<k-1> - p(a<k> - a<k+1>)²

	//calcula pi
	mpf_add(pi, sh_mem->a2, sh_mem->b2);
	mpf_mul(pi,pi,pi);
	mpf_div(pi,pi,sh_mem->t);
	mpf_div_ui(pi,pi,4);	

	f=fopen("saida.txt","w");
	gmp_fprintf (f, "%F1.10000000f\n",pi);

	mpf_clear(pi);
	destroy_shared_mem(sh_mem);

	return 0;
}


